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ABSTRACT 

An efficient algorithm has been devised to compute the crossovers in satellite altimetry. The significance of 
the crossovers is twofold. First, they are needed to perform the crossover adjustment to remove the orbit 
error.Secondly, they yield important insight to the oceanic variability. Nevertheless there is no published algorithm 
to make this very time consuming task easier, which is the goal of this note. The success of the algorithm is predi- 
cated on the ability to predict (by analytical means) the crossover coordinates to within 6 km and 1 second of the 
true values. Hence, only one interpolation/extrapolation step on the data is needed to derive the crossover coordi- 
nates in contrast to the many interpolation/extrapolation operations that are usually needed to arrive at the same 
accuracy level if deprived of this information. 

1. INTRODUCTION 

In satellite altimetry (in which the sea level relative to a reference ellipsoid that best approximates the shape 
of the earth is measured along the satellite ground track, e.g., Wunsch and Gaposchkin, 1980), the term "crossover" 
refers to the intersection of two ground tracks. The coordinates of a crossover point are comprised of the location 
(i.e., latitude and longitude) and the two times when the satellite passes over the crossover point. The crossover 
difference is the difference between the two sea level measurements. The crossovers are important in two aspects. 

First, the crossover difference reveals the temporal variability of the sea level and also avoids the geoid uncer- 
tainty (the geoid is an equipotential surface of the earth gravity field, to which a motionless ocean would conform). 
The sea level varies with respect to the reference ellipsoid primarily because of the earth gravity field and secon- 
darily due to the dynamical effects of the ocean. Presently the uncertainty of the geoid is at least comparable if not 
larger then the effects of ocean dynamics (e.g., Tai, 1983). Analyses based on crossover differences have shown 
and are continuing to reveal valuable insight on the sea level temporal variability (e.g., Cheney and Marsh, 1981; Fu 
and Chelton, 1985; Cheney et al., 1986.) 

Secondly, the crossovers play a pivotal role in removing the satellite orbit error (i.e., the uncertainty of the 
satellite's altitude) by the so-called crossover adjustment method (e.g., Tai and Fu, 1986; Tai, 1987). The orbit 

i 

error, which usually causes the crossover difference to be over 1 m in size, cannot be takeii lightly lest the crossover 

difference should reveal not the sea level temporal variability but the variability produced by the orbit error. Furth- 
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ermore, the basin-scale geoid uncertainty is small enough for the basin-scale circulation to be investigated (e.g., Tai 
and Wunsch 1983, 1984) and the orbit error reduction is of paramont importance in this case. 

For the exact repeat orbit, one can use along-track differences instead of crossover differences to reduce the 
orbit error if one is only interested in the temporal variability. However, the crossover adjustment has to be per- 
formed if one is interested in the absolute topography because ascending and descending tracks may have different 
orbit error characteristics which are only evident in crossover differences. 

There are many error sources in satellite altimetry (e.g., Tapley et al„ 1982). Thus, it is imperative that any 
estimate be deduced from many samples so as to minimize the effect of a particular error source or a particularly 
bad error realization because the errors are more or less independent of each other; and it is not unusual to find 
thousands or even millions of crossovers being treated in a single problem (e.g.. Marsh et al., 1982; Rapp, 1983). 
However, the simple prelude of finding the crossovers can consume more computing resources than solving the 
problem itself if an inefficient algorithm is used. Realizing the importance of a good algorithm, yet one would 
search in vain in the literature for a published algorithm to compute the crossovers. Since the launch of Geosat 
(March, 1985), the matter has taken on added urgency. The initial 18 months of the Geosat Mission is classified, but 
the crossover differences in this period are unclassified and would be available to the public if it were not such a 
cumbersome job. 

The purpose of this note is to make one good and tested algorithm available to the public. Hopefully a better 
algorithm will be published as the result of this note. In the following, it will be demonstrated that one can derive 
the coordinates of a crossover point to within 1 second and 6 km of the true values (1 km if an ad hoc formula is 
used) from the circular orbit approximation while compensating for the earth’s oblateness. Thus only one 
interpolation/extrapolation step is needed to derive the coordinates of the crossover from the data as oppose to three 
or four interpolation/extrapolation steps that are usually needed if deprived of this information. 

2. CIRCULAR ORBIT APPROXIMATION 

The orbit of the satellite is better described by an ellipse. However for the purpose of satellite altimetry, the 

i 

i 

orbit is made so circular that the circular approximation would greatly simplify the problem while incur little error. 

For example, the first eccentricity (see definitions in Section 3) of the Seasat’s orbit is merely 0.001, while its coun- 
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terpart for the earth is 0.0818 (i.e., one incurs much bigger error by assuming the earth is spherical. The effect of 
the earth’s oblateness will be discussed in Section 3). 

If one assumes the orbit is circular and the earth is spherical, the location of the ground track can be easily 
derived as a function of time. This has been done in Tai and Fu (1986, see their eq. (3), Appendix 1, and Figure 1). 
We can adapt their formulation to the present case; and the relevant equations are 

<j> = ./‘ sin -1 (rinAtsini) , (1) 

AX = tan -1 (tanAlcosi) - ClAt , (2) 

where the time, t, is nondimensionalized to make 2k correspond to the duration of one revolution; <|> is the latitude 

(the overhead bar conveys the fact that it is the geocentric latitude, see Figure 1); X is the longitude and defined to be 
in the range of [0, 2n); i is the inclination angle; £1 is the earth rotation rate relative to the orbital plane. Quantities 
with a subscript o are related to the equator crossing, and At = t-t„ , AX = X - X» (i.e., f„and X„ are the equator cross- 
ing time and longitude respectively). The adaptation is done so that equations (1) and (2) are valid for either an 
ascending or a descending track with j=l if ascending and j=-l if descending. Thus, the valid range of each vari- 
able is: -i £ <j> £ i, )At) <, kJ2, and |AA| 5 (l+k£l) it/2, where k=l if i > Jt/2 (i.e., retrograde), and k= -1 if i < tc/ 2 (i.e., 
prograde). 

2a. Determine the equator crossing time and longitude 

As a first step, the algorithm requires that the data be sorted into ascending and descending tracks, and the 
equator crossing time and longitude be determined for each track. Because data gaps (e.g., over land) often prevent 
the direct determination of these coordinates from data, there is a need to determine them analytically. From a point 
(with cordinates t, <j>, X) along the track, one can determine from (1) that for the equator crossing time, 

At = j sin~ l (sin fy/sini) , (3) 

t 0 =t-At , (4) 

and for the equator crossing longitude, from (2) and (3), 

X„ = \-tan~ l (tan A/cosi) + QAf . (5) 

Note that if A* should lie outside [0, 2 it), one can add or substract 2n from X 0 to make it fall in this range. 

2b. Determine the coordinates of a crossover point '' 

(1). Longitude 

\ 


Let a subscript a ( or d) convey the meaning of ascending (or descending). Then from the geometrical sym- 
metry entailed in equations (1) and (2), it is easy to see that the longitude of the crossover point must lie right in the 
middle between the two equator crossing longitudes of the two tracks. To be more specific. 


i X oa + X 0 j 

2 

* |X M -xj<;(i-o)7i . 

(6a) 

* ^ oa \ od , 

X = ±71 

2 

if |X M - X^l £(1+Q)n , 

(6b) 


where the sign in (6b) is such that (KX < 2 jc. And if (l-Cl)n: < |X M -X OI< | < (l+£2)rc, we have two possibilities: 

(1) two crossover points if i > tc/ 2, i.e., 

X, = bsL^SL , X a = X,±n (6c) 

(ii) no crossover point if i < n/2. 

(2) Time 

Knowing X (therefore AX), one can solve for At using eq. (2), which can be transformed to a more convenient 
form that avoids the evaluation of tan -1 , i.e., 

/(A/) = tan(AX + fiA/)- fanAfcosi =0 . (7) 

Hence, it becomes a problem of finding the zero of the transcendental function f. One can use either the ascending 

or the descending equator crossing longitude to form AX. Care must be taken (i.e., add or substract 2n) to make sure 

that |AX| <, (l+k£2)rc/2. Also note that if AX < 0 (or AX > 0), At lies between 0 and kn/2 (or -kn/2). Remember that 

k=l if i > n/2 and k=-l if i<n/2. 

(3) Latitude 

The geocentric latitude can be derived immediately from eq. (1). 

\ 

3. EFFECTS OF THE EARTH’S OBLATENESS 

As mentioned in the previous section, the earth is far more elliptical than the orbit. The effects are twofold 
(see Fig. 1). First, the subsatellite point (marked as g in Fig. 1, where the line joining it with the satellite is perpen- 
dicular to the local earth surface) is different from the point (marked s in Fig. 1) where the line joining the satellite 
with the earth center intersects the earth surface. Secondly, the coordinate data are given in terms of the geographic 
latitude (defined as the angle between the local vertical and the equatorial plane, i.e., the latitude that can be deter- 
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mined from local astronomical observations), which is to be distinguished from the geocentric latitude because of 
the earth’s oblateness (see Fig. 1, varibles with an overhead bar are geocentric, while those without it are geo- 
graphic). The subscript g (or s) distinguishes quantities related to point g (or s). 

The relation between <)> and <j> can be derived from (see Appendix) the following equation. 


jan<t> = P =n _ _ 1 
tan<|> U ’ \+e ' 2 


( 8 ) 


where a and b are iespectiely the major and minor axis of the earth; e is the first eccentricity [defined as 
e=(l -b 2 la 2 ) U2 ]\ e' is the second eccentricity [defined as e'=(a 2 /b 2 -l) 112 ]. These values are given in Table 1. 


The derivations in Section 2 all refer to <J>, , but the data refer to 4>, . Their relationship can be derived by solv- 
ing the following transcendental equation (see Appendix), 

, _ (9) 

ccos^ -bcosP b 

where c is the distance of the satellite to the center of the earth and p is related to [see eq. (10) below]. Given p, 
<j>, can be solved, or vice versa; and the relation between p and ij) f is given by 


tan<{> r = ^ tan ft . (10) 

Two inequalities are very helpful in specifying the range that the solution of equation (9) must fall in. First, to 
determine the equator crossing time and longitude, <j>, is given, but <j>, is needed in equations (3) and (4). Thus, to 
determine ^ from <J» f , <(>, is in the range (see Fig. 1), 

<i>, <$,<<!>, • (ii) 

Secondly, to determine the latitude of a crossover point, equation (1) yields the value of <[>,, while <|>, is the final 
answer. Hence, to determine <|>, from , (j> t is in the range (see Fig. 1), 

<{>/<<)>,< <t>, • (12) 

It has been discovered that there is no need to solve equations (9) and (10). The following two approxima- 
tions, representing the midpoints of the ranges specified by equations (11) and (12), actuall yield better results than 
equations (9) and (10). 


= OS ^ "(■ )/2 , 

$1 “(< 1 >j +< l >#)/2 . 


(13) 

(14) 
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The better accuracy (1 km versus 6 km, see Section S) does not portent physical consistency. Actually, the 6 

km accuracy in location is consistent with 1 second accuracy in time (i.e„ the ad-hoc approximation is simply a 

lucky guess). The orbit eccentricity (with can cause deviations from the circular orbit by typically a few parts per 

SMlrticVv 

thousand) is perhaps the cause here. 

To be physically more consistent, one could take the orbit eccentricity and the movement of the perigee (i.e., 
the point on the orbit ellipse which is closest to the center of the earth) into account But this would defeat the pur- 
pose of an efficient algorithm because one only needs to get close enough so that one interpolation/extrapolation 
step would give us the precise values. 

In summary, the oblateness only affects the latitude. Hence, it has no effect on the determination of the times 
and longitude of a crossover point (Section 2b) except that it may affect the equator crossing time and longitude if 
these coordinates have to be derived from the latitude (Section 2a) in lieu of the data. However, the oblateness 
effects can be accounted for by equations (9) and (10), or equations (13) and (14). As such, an accuracy to within 1 
km can be achieved. 

4. ALGORITHM 

A. Sort data into ascending and descending tracks. 

B. Determine the equator crossing time and longitude for each track. 

It is recommended that each category of data be assembled into an array, e.g., all the latitudes of a track can 
form a latitude array. Then there is the so-called table look-up (look-down) subroutine, which searches a strictly 
increasing (decreasing) array for the location index that a specified level is penetrated for the first time. Thus, speci- 
fying zero for the latitude array, one can quickly find the closest points to the equator to facilitate the 
interpolation/extrapolation to the equator. However, if the nearest points are too far away, one could apply equa- 
tions (3), (4), (3), and (13) to the nearest point to find the time and longitude of equator crossing. 

C. From two tracks, compute X (longitude of the crossover point) from eq. (6). ^ 


\ 
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Note that \ is the perfect candidate for the specified level in the table look-up (look-down) subroutine in step 
E except that the longitude array is not strictly increasing (descreasing) because of the jump at zero longitude. 

D. Compute At from eq. (7). Then compute /„and t d . 

E. Use r„and t d as the specified levels in the table look-up to find the nearest points in the data to the crossover 
point. 

F. Interpolation/extrapolation step. 

However, if the nearest points are not close enough, this step should be abandoned because the crossover 
difference would be contaminated by the geoid. 

5. RESULTS AND DISCUSSIONS 

The algorithm outlined in the previous section has been tested using the Geosat exact repeat data. The first 
descending track of the second repeat period (17 days) is selected as the basis. Then crossovers between this and 
the subsequent 244 ascending tracks (i.e., to the end of the 17-day period) are computed using the algorithm. The 
results of the first day (actually the first file. There are 17 files for each repeat period) are tabulated in Table 2. Note 
that 260 crossovers are generated in 17 days because 16 ascending tracks form two crossovers each with the 
selected descending track (see Appendix 2 of Tai and Fu, 1986, for explanation). It would be rather cumbersome to 
list them all. In Table 3, statistics of the differences between predicted (i.e., analytic approximaton) and true (i.e., 
deducted from data) values are presented. Note that the data only yield 147 crossovers (out of a theoretically possi- 
ble 260) because of data gaps. 

The descending track descends from (72.054°N, 284.959°E) to (72.054°S, 92.9 10°E) and crosses the equator 
at 189.134°E on 59877839.064 seconds counting from the beginning of the Geosat Mission (hereafter 59 million 
seconds are substracted to make the presentation easier). In Table 2 and 3, all unprimed variables are derived from 
the data. Primed variables are derived from the equations. 4*' is determined from equation (14) (i.e., a simple aver- 
age), while 4>" is derived from equations (9) and (10). As discussed in Section 3, <}>' approximates <f> better than <(>". 
From Table 2, <]>' is accurate to within 0.01° and <(»" is accurate to within 0.06°. The longitude is highly accurate 
except near high latitudes where the tracks are going more east-west than north-south. But even there, th& 
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maximum error is only about 0.01°. The time discrepancies are generally less than 1 second except for two cases 
where they are around l.S seconds. Note that the orbit eccentricity is partially accounted for here by adapting an 
orbit period of 6031.4 seconds north of the equator, while using 6043.6 seconds south of the equator. Because the 
purpose is to compute the coordinates of the crossover point, all available data have been used (i.e., the land values 
are included). These general remarks are substantiated by Table 3. 

To compute all 260 crossovers, it takes about 0.5 second cpu time of the San Diego Supercomputer Center’s 
Cray XMP-48. It translates to about two-minute cpu time to compute all the crossovers in the 17-day period. 
Because the predicted values are so close to the real values, only one interpolation/extrapolation step is needed 
versus the situation that many steps are needed if starting from scratch. Without any optimization, the prediction 
step takes less than one third the time that is needed for the search and interpolation/extrapolation step. 

The transcendental equation in the prediction step is solved by a method, which combines the bisection and 
the secant rule. One could easily optimize the prediction step by solving and storing the solutions to be used later as 
the starting points for the interations. For instance, eq. (7) can be solved for At with AX ranging from 1° to 192° 
with 1® increment. Then when one gets a AX of 99.5°, he can use the solutions for 99° and 100° to start the intera- 
tion. Furthermore, because the location can be predicted so accurately, when the crossover point is obviously over 
land, the extremely time consuming search step (i.e., table look-up, look-down) can be avoided altogether, for 
example, by setting up 1° by 1° land flags. 


APPENDIX 


To get the relation between $ and $, an ellipse can be expressed as 


2 v 2 

— 4 - — = 1 

2 b 2 ' 


or as 


Thus from (15), (16), 


x = a cosp , 
y = b sinp . 



tanp 


(15) 

(16) 


(17) 

\ 


and 


i.e., the tangential unit vector at the point (x,y) and in the direction of increasing p is proportional to (-a sinp, b 
cosP). Thus the normal unit vector is proportional to (b cosp, a sinP). Hence, 

tan<j> = -^-tanp (20) 

b 

From (17) and (20), eq. (8) is the result To derive eq. (9), let us consider Fig. 1. The coordinate of the satellite is (c 
cos <j>, , c sin (j>,). The coordinate of point g can be expressed as (a cosp, b sinP) [see equations (15), (16)]. thus, one 
can form the vector from point g to the satellite and eq. (9) would result if one considers the ratio of the y com- 
ponent to the x component of this vector. 
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FIGURE LEGEND 


Fig. 1. 


An exaggerated view of the earth and the satellite. 
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Table 1. The relevant Geodetic Reference System 1980 values adopted in this paper and certain orbit parame- 
ters 


a = 6378137 m 


semimajor axis 


b = 6356752 m 
e 2 = 0.0066943800229 
e * = 0.00673949677548 
c = 7163 km 


semiminor axis 
e = first eccentricity 
e' = second essentricity 

distance between the satellite and the center of the earth 
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Table 2. Computed and predicted crossover coordinates where the time is in seconds, latitude and longitude in 
degrees. 


X 0a 

td 

t'd 

‘a 

t 'a 

X 

X' 

♦ 

<t>' 

♦" 

356.581 

N.A. 

876395.66 

N.A. 

882304.23 

N.A. 

272.858 

N.A. 

71.664 

71.709 

331.499 

876465.92 

876465.69 

888271.47 

888271.69 

260.319 

260.317 

70.406 

70.406 

70.453 

306.417 

N.A. 

876548.71 

N.A. 

894225.99 

N.A. 

247.776 

N.A. 

68.000 

68.052 

281.335 

N.A. 

876657.10 

N.A. 

900155.21 

N.A. 

235.235 

N.A. 

63.813 

63.872 

256.254 

876815.08 

876814.40 

906034.89 

906035.60 

222.697 

222.694 

56.512 

56.506 

56.575 

231.172 

877067.34 

877066.59 

911820.25 

911821.03 

210.155 

210.153 

43.370 

43.360 

43.435 

206.090 

877477.72 

877477.28 

917447.45 

917447.91 

197.612 

197.612 

20.560 

20.551 

20.601 

181.007 

878017.13 

878017.66 

922945.60 

922945.05 

185.071 

185.071 

-10.149 

-10.146 

-10.172 

155.925 

878487.38 

878488.84 

928512.87 

928511.37 

172.529 

172.530 

-36.588 

-36.584 

-36.656 

130.843 

878789.53 

878790.95 

934248.23 

934246.78 

159.988 

159.989 

-52.739 

-52.740 

-52.812 

105.761 

878974.18 

878975.28 

940101.13 

940100.01 

147.444 

147.448 

-61.698 

-61.700 

-61.762 

80.679 

N.A. 

879097.57 

N.A. 

946015.34 

N.A. 

134.907 

N.A. 

-66.778 

-66.832 

55.597 

879187.41 

879187.90 

851963.24 

951962.64 

122.354 

122.366 

-69.712 

-69.713 

-69.762 

30.515 

879261.35 

879261.63 

957926.87 

957926.58 

109.813 

109.825 

-71.328 

-71.330 

-71.375 

5.434 

879327.39 

879327.47 

963898.32 

963898.17 

97.274 

97.284 

-72.006 

-72.008 

-72.052 

5.434 

N.A. 

876372.60 

N.A. 

966853.04 

N.A. 

277.284 

N.A. 

71.894 

71.938 
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Statistics of the differences between predicted (primed) and real (unprimed) values. Time is in 
seconds. Latitudes and Longitudes are in degrees. 


Statistics 

t 'a~ta 

t 'i-t d 

r-x 



Maximum 

0.79 

1.53 

0.011 

0.0058 

0.0667 

Minimum 

-1.56 

-0.78 

0.078 

-0.0101 

-0.0746 

Mean 

-0.34 

0.34 

-0.0004 

-0.0023 

-0.0162 

Standard 

Deviation 

0.82 

0.79 

0.0098 

0.0041 

0.0532 

RMS Value 

0.88 

0.86 

0.0098 

0.0047 

0.0555 




